function [betas_tightness, betas_f] = reg_hpxr_real(gesol,Phi_NLz,prob_Nxz,FLAG_RUN_PARALLEL)
% slope coefficients for the following regression:
% hpxr(n,t+12) = a + b*X(t) + e(t)
% hpxr(n,t+12) = p(n-1,t+12) - p(n,t) + p(1,t)
% X(t) is either tightness(t) or f(t)
% maturity is n=2,3,4,5 years

    if nargin<4; FLAG_RUN_PARALLEL = false; end

    mean_tightness = sum(prob_Nxz(:).*gesol.theta(:));
    var_tightness = sum(prob_Nxz(:).*(gesol.theta(:).^2)) - mean_tightness^2;
    mean_f = sum(prob_Nxz(:).*gesol.f(:));
    var_f = sum(prob_Nxz(:).*(gesol.f(:).^2)) - mean_f^2;

    % only compute matrix powers once (this is very slow)
    Phi_NHz_12mo = mpower2(Phi_NLz,12);

    betas_tightness = zeros(1,4);
    betas_f = zeros(1,4);

    idx_1y = find(gesol.real_bonds.maturities==12);

    [NN, NL, Nz] = size(gesol.theta);
    temp_NLz = zeros(NN*NL*Nz,1); % for storage

    temp_NLz(:) = prob_Nxz(:).*reshape(log(gesol.real_bonds.prices(:,:,:,idx_1y)),[],1);
    mean_p1 = sum(temp_NLz);
    mean_p1_tightness = sum(gesol.theta(:).*temp_NLz);
    mean_p1_f = sum(gesol.f(:).*temp_NLz);

    if FLAG_RUN_PARALLEL
        
        theta_flat = gesol.theta(:);
        f_flat = gesol.f(:);
        prob_Nxz = prob_Nxz(:);
        p_buy_list = cell(1,4);
        p_sell_list = cell(1,4);
        
        for n = 2:5
            
            idx_pnbuy = find(gesol.real_bonds.maturities==12*n);
            idx_pnsell = find(gesol.real_bonds.maturities==12*(n-1));
            
            p_buy_list{n-1} = gesol.real_bonds.prices(:,:,:,idx_pnbuy);
            p_sell_list{n-1} = gesol.real_bonds.prices(:,:,:,idx_pnsell);
            
        end
        
        parfor n = 2:5

            temp_NLz = prob_Nxz.*reshape(log(p_buy_list{n-1}),[],1);
            mean_pnbuy = sum(temp_NLz);
            mean_pnbuy_tightness = sum(theta_flat.*temp_NLz);
            mean_pnbuy_f = sum(f_flat.*temp_NLz);

            temp_NLz = prob_Nxz.*(Phi_NHz_12mo*reshape(log(p_sell_list{n-1}),[],1));
            mean_pnsell = sum(temp_NLz);
            mean_pnsell_tightness = sum(theta_flat.*temp_NLz);
            mean_pnsell_f = sum(f_flat.*temp_NLz);

            mean_hpxr_n = mean_pnsell - mean_pnbuy + mean_p1;
            mean_tightness_hpxr_n = mean_pnsell_tightness - mean_pnbuy_tightness ...
                + mean_p1_tightness;
            mean_f_hpxr_n = mean_pnsell_f - mean_pnbuy_f ...
                + mean_p1_f;

            betas_tightness(n-1) = (mean_tightness_hpxr_n - mean_tightness*mean_hpxr_n)/var_tightness;
            betas_f(n-1) = (mean_f_hpxr_n - mean_f*mean_hpxr_n)/var_f;

        end
        
    else
    
        for n = 2:5

            idx_pnbuy = find(gesol.real_bonds.maturities==12*n);
            idx_pnsell = find(gesol.real_bonds.maturities==12*(n-1));

            temp_NLz(:) = prob_Nxz(:).*reshape(log(gesol.real_bonds.prices(:,:,:,idx_pnbuy)),[],1);
            mean_pnbuy = sum(temp_NLz);
            mean_pnbuy_tightness = sum(gesol.theta(:).*temp_NLz);
            mean_pnbuy_f = sum(gesol.f(:).*temp_NLz);

            temp_NLz(:) = prob_Nxz(:).*(Phi_NHz_12mo*reshape(log(gesol.real_bonds.prices(:,:,:,idx_pnsell)),[],1));
            mean_pnsell = sum(temp_NLz);
            mean_pnsell_tightness = sum(gesol.theta(:).*temp_NLz);
            mean_pnsell_f = sum(gesol.f(:).*temp_NLz);

            mean_hpxr_n = mean_pnsell - mean_pnbuy + mean_p1;
            mean_tightness_hpxr_n = mean_pnsell_tightness - mean_pnbuy_tightness ...
                + mean_p1_tightness;
            mean_f_hpxr_n = mean_pnsell_f - mean_pnbuy_f ...
                + mean_p1_f;

            betas_tightness(n-1) = (mean_tightness_hpxr_n - mean_tightness*mean_hpxr_n)/var_tightness;
            betas_f(n-1) = (mean_f_hpxr_n - mean_f*mean_hpxr_n)/var_f;

        end
    
    end

    % multiply by 100 since hpxr is not in percentages
    betas_tightness = 100*betas_tightness;
    betas_f = 100*betas_f;

end